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I. INTRODUCTION 

There has been a lot of interest in the Kitaev model 
in recent years. The quantum mechanical spin-1/2 Ki- 
taev model is exactly solved in 2D. The Hamiltonian can 
be diagonalized exactly in terms of Majorana fermions 
The model exhibits a phase transition from a phase 
with finite correlation length to one with long-range cor- 
relations as the ratios of coupling constants in different 
directions are varied 0,111. It has topological excitations, 
and their robustness with respect to noise makes it an in- 
teresting candidate for quantum computing 

In a recent very interesting paper, Baskaran et al. 
studied a generalization of this model with spin-5' at each 
site and identified mutually commuting Z2 variables that 
are constants of motion for arbitrary S For large S, 
the spins can be approximated as classical 0(3) vector 
spins. Baskaran et al. showed that the classical ground 
state of the model has a large degeneracy. They argued 
that though a naive averaging over these ground state 
configurations would suggest that the system is disor- 
dered at zero temperature, for large S, the quantum 
fluctuations of spins have lower energy for a subset of 
the classical ground states. These states get more weight 
in the quantum mechanical ground state, and the quan- 
tum model shows long-range order in the ground state, 
an example of quantum order-by-disorder. 

The finite-temperature fluctuations in the classical 
model behave qualitatively like the zero-point fluctua- 
tions in the quantum model, and it is interesting to ask 
if temperature fluctuations can induce order-by-disorder 
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in the classical Heisenberg spins with Kitaev couplings, 
just as the quantum fluctuations are expected to, in the 
large- quantum model studied by Baskaran et al.. In 
this paper we point out that there is a qualitative dif- 
ference between the classical and quantum mechanisms 
of order-by-disorder. For the classical Kitaev model, the 
contribution of nearby states to the restricted partition 
function in the limit of very small temperature, with 
states summed only over the neighbourhood of a given 
classical ground state, is exactly the same for almost all 
ground states. For the ensemble of ground states, we 
establish an exact equivalence to a solid-on-solid model 
with nearest neighbour coupling at a finite temperature. 
We argue that bond-energy bond-energy correlations de- 
cay as R"^ with distance R at zero temperature. Our 
Monte Carlo simulations also support this conclusion. 
We derive upper and lower bounds on the ground state 
energy of the quantum spin-S" model. We also study the 
ground state energy of the quantum system using varia- 
tional wavefunctions, avoiding the divergences present in 
the spin- wave expansion. 

Kitaev's pioneering work has led to a large amount of 
further research. The spectrum of the different phases 
of this model have been extensively studied [6]. Propos- 
als for experimentally realising this model using polar 
molecules and ultracold atoms trapped in optical lattices 
have recently been made 0, Q- Alternate methods of 
solving this model usiiig Jordan- Wigner transformations 
have been proposed while perturbative studies have 
also proved fruitful [10|, |n| . The Kitaev model has also 
been studied on other two dimensional lattices 
Exact solutions have been obtained for the Kitaev model 
on certain three dimensional lattices [3 E^- There is 
also a fair amount of earlier work on order-by-disorder 
in classical systems [2^ [2l[. It has been studied a lot 
in the context of magnetic systems with frustration, such 
as spin systems with nearest nei ghb our antiferromagnetic 
interactions on different lattices [2^ [l^ . 
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II. DEFINITION OF THE MODEL 

We consider classical Heisenberg spins on a hexagonal 
lattice. We consider a finite lattice, with periodic bound- 
ary conditions. There are L hexagons in each row, and 
M rows of hexagons [L and M both assumed even). The 
total number of hexagons is iV = LM, and the number 
of sites is 2N . The bonds of the lattice are divided into 
three classes, X, Y and Z, according to their orientation 
(Fig. [T]). The hexagonal lattice consists of two sublattices 
denoted by A and B. We label sites in the A sublattice 
by a{l,m) and the corresponding B sublattice site con- 
nected to it via a Z bond by b{l,m). We define three 
bond vectors Cx-, Ry, &z as the vectors from any A site to 
its three neighbours via the X, Y and Z bonds respec- 
tively (Fig. [1]). Thus we have a(Z, m) + = b{l, m — 1), 
a{l, m)+ey = b{l+l, m—1) and a{l, ra)+ez — b{l, m). We 
define the hexagonal plaquette {I, m) to be the hexagon 
whose topmost point is a{l, m). A bond will be specified 
by the {I, m) coordinate of its A-lattice end point and its 
class X, Y or Z. For instance, (a(Z, m); x) = {l,m]x) 
is an X-bond with a{l,'m) at one of its ends. The pe- 
riodic boundary conditions are implemented by making 
a{l,m) = a{l + L,m) ^ a{l ~ ^,m + M). 

At each lattice site i there is a three dimensional vector 
spin S^ = {S{',Se', S^^) of unit magnitude. Thus S'^^^ + 
Si^"^ + Si^'^ = 1 at every site. The Hamiltonian of the 
system is given by 



ya(l-l.m+l) j^a(l,m+l) 
^b(l-l,m) X^(l,m) l)'b(l+J,m) 



a(l+],m) 




FIG. 1: (Color online) A hexagonal lattice depicting the la- 
belling scheme for sites, and the x, y and z bond classes. Sites 
in the A- and B- sublattices are denoted by filled and open 
circles respectively. 



evaluated as 



exp[-/3 cos + 5,^(,_i,„.+i)' + ^,^(,,^)' 



(3) 
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cy qy 
'-'a'-'a+Cy 



^a^a+ej (1) 



where 9 is the angle between the vector b{i,m) and i^. 
This immediately yields 



As the hexagonal lattice is bipartite, for classical spins, 
without loss of generality we take J > 0. The Hamilto- 
nian does not have rotational symmetry in the spin space, 
but it has a local symmetry : for any bond {l,m;a), 
the Hamiltonian is invariant under the transformation 



a.{l,' 



ca qa 



_ qa 



Now 



sinhf/3J5^, ^,.'^ + 8^,, , + ,^1 

y a(t,m+l) a(t— l,m+l) a{l,m) J 

qx 2 I qy A_ Qz 2 



n %^ n 



(4) 



(5) 



(i,m) 
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III. FINITE TEMPERATURE PARTITION 
FUNCTION 



We thus obtain the effective Hamiltonian for the spins 
on the A-sublattice alone as 



The partition function of the system at finite temper- 
ature Z[f3] is given by 

^[/^] = / n(^) cxph/3if] (2) 

where = T. The index s runs over all sites of 

the lattice. The integral over each B-site is of the 

form Wi^ra = J cf^ b{i,7n) exp[-i3~i b{i,7n)-~^] where ~f' = 
Sa(Urn+i)i + Sa{i-i^m+i)3 + Sa(L,n)k. This Can in turn be 



H,ff{{Sa},l3) = 

~ o X! ^(/5(\/'5'a(i,m+l)^ + ^l(l-l,r,i+l)^ + '5'a(i,m) ^) ) ' 
(i,m) 

(6) 

where 

^, , , .sinh(a;), 

F{x) = log[— (7) 

For a given configuration of spins {S"}, to each bond 
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(Z,m;a) of the lattice, we assign a vector e(l,m;a)ea, 
with e{l, m; a) given by 



(8) 



We define the discrete divergence of the e-field on the 
A, and B sublattices as 

V.e = e(a(^ m); a) at site a{l, m), (9) 
V.e = ^ £{b{l, m) — Ba, a) at site h{l, m). (10) 



Clearly, the divergence of the field e at any site on the 
A-sublattice is 0. In addition, for ground state config- 
urations ( the proof of this assertion is given the next 
section), even for bonds meeting at any site, G i?, we 
have 



e(6(/, m) — e^; a) = 0, for all sites b ^ B. (11) 



For non-ground state configurations, the sum of e vari- 
ables at each i?-site is no longer unity. We parametrise 
this deviation using the new variable Q defined by 



ea;a). 



(12) 



We can think of these variables Qb(i,m) as charges 
which are placed on the sites of the B-sublattice; there 
are no charges on the A-sites. We define Qa{i,m) = 0, for 
all sites a(/, m). Given the values of these charges, we can 
construct the corresponding electrostatic potential field 
(h defined at all sites s of the lattice such that 



-Qs\ for all sites s. 



(13) 



Here is the discrete laplacian on the lattice. These 
equations can be solved explicitly, and determine the po- 
tential (/)(s) completely, up to an overall additive con- 
stant, so long as the total charge in the system is zero. 
Explicitly, we have 



Y.G{s,s')Q, 



(14) 



where G{s, s') is the lattice Green's function. Then, as 
\/^(f>{s) = —Qs — — V.e(s), we see that e -I- Vcfi has no 
divergence and can be expressed in terms of the curl of 
a new field /. We define the scalar field /(/, m) = /;.,„ 
attached to the hexagons of the lattice (as shown in Fig. 
[5]) such that the difl[erence in the / field between two 
neighbouring plaquettes is equal to the value of e -I- V0 
along the shared bond. This satisfies the divergence-free 
condition for the field e -|- Vcj). Let s be any site on the 
A-sublattice. Its neighbours are sites s + Bx-, s + By, s + 



Cz . Let the three hexagons to which s belongs be hi , /12 
and /13 (Fig. [2]). If the site s = a{l,m), then hi will 
have coordinates (Z — 1,to + 1), and similarly for other 
hexagons. Then, for all sites s, the /- field is defined by 



e{s,x) = (t>{s) - i 
e{s,z) = (p{s) - 



{s + ex) + f{hi)-f{h2) 
{s + ey) + f{h2)-f{h3) 
{s + e,) + fih3)~fihi) 



(15) 



Given the fields e{s,a) and (/i(s), we assign any fixed 
value to /(Z, m) at one particular hexagon, then the value 
of the /-field at neighbouring hexagons is completely de- 
termined. Thus for a given configuration of {5"}, we 
can determine the /-fields at all hexagons up to an over- 
all additive constant. 



O' 



s+e 




s+e^ 



FIG. 2: Figure depicting the definition of the e and h variables 
on the bonds and plaquettes around an A-site s. 

We will use the values of {Qs} and {/(/,m)}, instead 
of {5"^} to specify the spin-configurations. The number 
of variables S""^ are 3iV in number, with N constraints 
between them, thus there are 2N independent real vari- 
ables. As the variables Qs satisfy the constraint J^s Qs — 
0, there are N — 1 independent parameters Qs- Also, 
there are only (N — 1) independent parameters f{l,m), 
as these are defined only up to an overall additive con- 
stant. We need two additional linearly independent vari- 
ables to complete our new set of coordinates. We choose 
these to be Ri = m; z) + (j){b{l, m)) — <j){a{l, m)) 

and i?2 = J2m *^(^' y) + 'f'iKl + 1, TO - 1)) - (/>(a(Z, m)), 
which lead to 



f{l + L,m) = f{l,m) + Ri 
M 

f{l~ —,m + M)^ f{l,m) + R2 



(16) 



Assuming that the </)- field, obtained in Ea.([T^ is pe- 
riodic on the torus, with (/){a{l,m)) = (j){a{l + L,m)) — 
(t){a{l — M/2, m + M)), Ri is independent of m (and R2 
of I), and these correspond to fixing the boundary condi- 
tions for the /(/, m). 
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The condition 2/3 > e{s,a) > —1/3 at each bond 
impUes constraints on the ahowed range of f{l,m) and 
Q{l,m). Now, given the values of the Q, f and (j) fields, 
one can systematically reconstruct the e field (and thus 
the spin configuration) . The value of e at any bond in the 
bulk can be evaluated from the f's at the neighbouring 
plaquettes and the (f) at the end of the bond as shown in 
Eq. The e's at the edges of the lattice are deter- 

mined by the values of / at the plaquettes next to the 
edge, which can be obtained from i?i (or R2) using Eq. 
(|16p . Thus, in the allowed range, the transformation from 
{Q, f, R} to {S°'^} is invertible. We can now express the 
partition function in terms of these new variables. 

Firstly we analyse the phase space factors in the parti- 
tion function as we change variables from {S} to {/, Q}. 
The phase space integral for each A-site spin is 



ds = I ds''dsyds'S{s''^ + 



S^^-l). (17) 



We change our integration variables from to (5^ )^. 
We have dS^ = d{S^)^ /{2y/{Sfy) and similarly for the 
y and z components. The f's and Q's are linear functions 
of the (5")^'s {a = x, y, z), hence the Jacobian matrix of 
transformation for this change of variables is a 2N x 2N 
constant matrix. Also, the determinant is non-zero as 
the transformation is invertible. The partition function 
at finite temperature, up to an unimportant constant, is 
thus given by 



ZIP] = (Const.) 



/ 



dRidRi 



n (i + E(bond) 



-1/2- 



exp 



^F(Vi + Qi,(i,m)) 



(18) 



where, Tltonds denotes the product over all bonds 
{l,m,a) of the lattice, with the e variables defined in 
Eq. ©, and F{x) is defined by Eq. ©. 



IV. CHARACTERISATION OF THE GROUND 
STATE MANIFOLD 

Baskaran et al. defined Cartesian states as states 
where each spin is aligned in the direction along a Carte- 
sian axis (x-, or y- or z-) Q. We can construct a Carte- 
sian ground state of H, by constructing a dimer covering 
of the hexagonal lattice. Spins at the ends of a dimer 
of type a {a = X,Y or Z), are aligned parallel to each 
other in the direction a (either both having S" — +1, 
both having S"" = —1). This state has an energy ~N J . 
Then corresponding to a dimer covering, there are 2^ 
Cartesian ground states. The number of dimer coverings 
of the hexagonal lattice increases as 1.38''^ [l^l, hence the 
number of Cartesian ground states increases as 2.76^. 
Baskaran et al. also showed that for any two Cartesian 



states, there is a one-parameter family of ground states 
that connects them, thus forming a network of ground 
states. 

However, they did not provide a proof that no states 
of lower energy can be formed, or study other possible 
ground states. In this section we characterise the entire 
set of ground states of this model. 

In the large /? limit, F(x) can be replaced by x. There- 
fore the ground state energy Eq of the system is given 
by 



Eq = - J Max 



1,171) 



(19) 



As ^/x is a convex function of x, for all real positive 
a;^, i = 1 to iV, we have 



N 



(20) 



With equality holding only when all the equal. 
Using this inequality in p^ . we get 



ErmniiSa}) > 'JN. 



(21) 



This result can also be arrived at by noting that 
F[/3y/x] is a convex function for all /3,x > 0. As the 
equality sign in Eq. ([T^ holds only when all the terms 
within the square root are equal, the necessary and suf- 
ficient condition for the ground state configuration is 



Qb = 0, for all sites b G B. 



(22) 



Since the Q -field, and hence also the (^-field are ex- 
actly zero everywhere in the ground states, the manifold 
is described only by the /-field. Correspondingly, the 
equations ([T5|) simplify to 



€{s,x)=fihi)~f{h2) 

e{s,y)^f{h2)-.f{h^) 
e{s,z)=,f{h3)-f{tH) 



(23) 



The set of states forms an 1 dimensional manifold, 
parametrized by the variables {/}, with the boundary 
conditions on these given by Ri and i?2- It is a convex 
set whose extremal points correspond to the Cartesian 
states studied by Baskaran et al.. 

We define the restricted partition functions for a fixed 
{f{l,m)} by integrating over the Q's 



(l,m) ■' 



exp[-/3//e//[{/,Q}] 



(24) 
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For large /?, the integrand in Eq. (p4)) is sharply peaked 
at Qs = 0. We can use the method of steepest descent to 
find the value of this integral (integrating over the — 1 
Q variables). We expand -ffe// in a power series in Q's 

H,ff = Eo + Y^Qs^ + ... (25) 

s 

The linear term in Q in the function H^f / vanishes since 
the '5s = 0. While the range of the Qs integrals de- 
pend on {/i^m}, for large /3, when the width of the peak is 
much smaller than the range of integration, and the peak 
is away from the end points of the range, each integra- 
tion to leading order is independent of {/i.m} and gives 
a factor C/3~^/^ where C is a constant. The restricted 
partition function Z [{/},/?] in the limit of very small 
temperature to leading order in /3, is /3^'-^^^-'/^Zo[{/}]. 
Where Zo[{/}] is given by 



V. EQUIVALENCE TO THE SOLID-ON-SOLID 
MODEL 

We note that may be interpreted as the partition 
function of a solid-on-solid (SOS) model, with a real 
height variable located at sites {l,m) of a triangu- 
lar lattice and interacting via an effective Hamiltonian 
Hsos- This Hamiltonian depends on the temperature 
T of the spin model. At a finite temperature Hsos is 
determined by integrating over the Q variables in the 
restricted partition in Eq. ([21)) . Hsos{T > 0) has 
some weak long range couplings. However at T = 
the Q field is identically zero, which leads to a purely 
nearest-neighbour, but non-quadratic coupling between 
the height variables given by 

Hsos{T = 0) = - [VUi,^-i - //-l,m) 

(l,m) 

+ V{fl^^ - fLm-l) + Vifl.i,m - fl,m)]. (30) 



Zo[{/}]= hm /3(^-i)/2z[{/},/3] 



Const. 



n [i + ^(bond)] 



-1/2 



.bonds 



(26) 



Thus for all fixed {f{l,m)}, the integration over fluctua- 
tions in {Q} produces the same temperature dependent 
weight factor, in the limit of large /3. For evaluating av- 
erages in the limit of low temperatures, we can ignore 
the Q-degrees of freedom, and set them equal to zero. 
Now, the zero-temperature partition function, i.e.- the 
partition function in the limit /3 — > oo dcflncd as 



Zo = hm /3(^-i)/2Z[/3] 

p—^oo 



dRidR2 
can be expressed as 



n / df{l,m) 

[1,771) 



Zo[{f}] (27) 



Zo 



where 



dRidR2 



n / '^fi^'^) 

{l,m) 



exp 



(28) 



V{x) = i log(i + x), for - 1/3 < .T < +2/3; 
— +00, otherwise. 



(29) 



and the sum over denotes the summation over all 
nearest neighbour hexagons i and j. 



We note that the Hamiltonian has a term log(^ -I- 
e(bond)), which diverges when e(bond) tends to —1/3. 
Thus the Cartesian states of Baskaran et al. have a large 
relative weight, which has a divergent density. However, 
this divergence is an integrable divergence, and the ac- 
tual measure of the Cartesian states in the ensemble of 
states at zero temperature is zero. 

The gauge symmetry of the model has a consequence 
that all correlation functions of the type (S'^^S'fj) with 
sites si and S2 not nearest neighbours are zero [27| . 
The simplest nontrivial correlation functions, for non- 
neighbour si and S2 are of the type ((>S'"i)^(5'f2)^). The 
convergence of the high temperature expansion of the 
partition function implies that these correlations fall ex- 
ponentially with distance, at small /?. As there is no 
phase transition as /3 oo, we expect this behaviour 
for all < /3 < oo as well, with the correlation length 
increasing as a function of /3. 

We now argue that, at zero temperature, this correla- 
tion function decays as for large separations R. 

The SOS model has the symmetry that changing all 
heights by the same constant leaves the Hamiltonian un- 
changed. Though the interaction is a strongly non-linear 
function of fb(i,m)+e^ - fb{i,m)+e^,, One expects that in 
the high-temperature phase of the SOS model, the long- 
wavelength hydrodynamical modes in the system will still 
be sound-like, with effective Hamiltonian |V/p, which 
gives rise to the spectrum given by oj^ cx fc^. For two 
sites si and S2 separated by a large distance R 



((/ai ^ log i?. 

This implies that 



(V/,,.V/.J 



1 



(31) 



(32) 
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Since the energy density variables (S'f )^ are propor- 
tional to V/ (at zero temperature), we conclude that the 
connected part of the bond-energy bond-energy correla- 
tion function 



1 



(33) 



At infinite temperature, the spins at different sites are 
completely uncorrelated. This is not true for the / vari- 
ables, which have non-trivial correlations even for /3 = 0. 
In the Appendix we calculate the leading behaviour of 
((/i?, — /o)^) at large i? for /3 = 0. We have, at infinite 
temperature 



2^/3 
457r 



Iog[i?] + O(l) for large R. (34) 



VI. MONTE-CARLO SIMULATIONS 

In this section we present results from Monte Carlo 
studies of this model for the zero temperature as well as 
for non-zero temperatures. 

We simulated the effective Hamiltonian H^f / (Eq. [6]) , 
obtained by integrating out spins on the i3-sublattice. 
For the finite temperature simulations, two kinds of 
moves were employed — single spin moves and cluster 
moves. 

We discuss single spin moves first. In any given state, 
we choose a site, a(Z, m). We generate a gaussian random 
vector r = {rx,ry,rz), whose variance is proportional to 
the temperature T. The proposed single spin move is 

to 



then to change the spin at site alljin) from 5*^ 
S'aii,m)^ given by 



{l,m] 



\Sa{l.r. 



(35) 



If the change in the effective Hamiltonian by the 
move is AiJ, the move is accepted with probability 
Min(l, e"''^^). Clearly, this satisfies the detailed bal- 
ance condition. 

While these single spin moves, in principle, are suffi- 
cient for correctly sampling the entire phase space, we 
also employed hexagon update moves to speed up the 
simulations at low temperatures. Given any configura- 
tion, we randomly choose a hexagon on the honeycomb 
lattice. To obtain the new configuration of spins we move 
along this hexagon, alternately adding and subtracting 
a quantity. A, to the bond-energies, and then comput- 
ing the spin components which give rise to these bond- 
energies. In Fig. [21 suppose the topmost A-site is si, 
then ei = S^^'^ — ^ and €2 = S^^'^ — ^. Now ei is changed 
to ei -I- A and 62 to £2 — A. This then fixes the new S^^ 
and 5*1^ (up to a randomly chosen sign), leaving S^_^ un- 
changed. Clearly, this leaves the sum of squares of the 




£2 -A 



e,-A 

o 





FIG. 3: The hexagon update move in the Monte Carlo sim- 
ulations. The e's depict the bond-energy variables associ- 
ated with each A-site (depicted by filled circles). A random 
number uniformly distributed between —a and +a is alter- 
nately added to and subtracted from the bond energies on 
the hexagon. The move is rejected if any of the bond energies 
fall outside the interval [— |, §]. 



spin components unchanged. This change is also made to 
the four other bonds on the hexagon [Fig. [3]. The value 
of A is chosen uniformly in the interval [—a, a], where a is 
a parameter. The proposed move is rejected if any of the 
bond-energies faU outside the interval [— 5, §]• Since the 
sum of bond-energies at each site (A and B) remains con- 
stant, these hexagon update moves leave the value of the 
effective Hamiltonian unchanged. These moves therefore 
play a crucial role in efficiently sampling the configura- 
tions close to the ground states at very low temperatures. 
We take the ratio of the phase space factors of the two 
states and accept or reject the proposed configuration ac- 
cording to the Metropolis rule. Clearly, this also satisfies 
the detailed balance condition. For the zero temperature 
simulations only the hexagonal updates were used. 

Monte Carlo simulation data presented in this section 
has been computed for L x L triangular lattices of A- 
sublattice spins of various sizes, with L ranging from 30 
to 256. 6 X 10^ Monte Carlo updates were made per site 
of which the first 6 x 10^ were not used in computing the 
correlation functions. Correlation functions were calcu- 
lated after every 6 updates per site. 

We calculated the correlation function C(f) = 
('^o(/ m)^^l(i' m')^) ~ \ ^o'" various lattice sizes and tem- 
peratures. Here r is the vector from site a(/, m) to 
a(r, m'). We find that this correlation decays quite fast, 
and at finite temperatures, is very small except for a few 
points around the origin. At zero temperature, C{r) is 
oscillatory along the (1,0) direction (and is periodically 
negative). We observe a clear behaviour along the 
direction, as plotted in Fig. El 

We looked for a signal of possible order in the ground 
state ensemble of the type proposed by Baskaran, et al.. 
One way of determining if there is any periodic order 
in the system is to study the structure factor which we 
define as the Fourier transform of C(r) 



(l',m') 



(36) 

where the summation is over all sites a(l' , m') for a fixed 
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r 



FIG. 4: (Color online) Plot of the zero temperature correla- 
tion function C(r) = (S'A^(O)SA^(r)) - | versus distance, r 
along the direction. These correlations follow a power law 
behaviour with exponent ~ —2. The line has a slope of —2. 



0.26 




r 



FIG. 6: (Color online) Graph of the zero temperature correla- 
tion function ((/(O) - f{r)f} versus distance, r, for different 
lattice sizes L, showing a log r dependence in accordance with 
the mapping to a height model at a finite temperature. 




where the r summation extends over all lattice sites. Two 
prominent peaks are visible at (— ^) and (^, — |^). How- 
ever, these peaks do not diverge with system size in our sim- 
ulations. 



a{l,m) with r as defined earlier. The structure factor 
S{k) would have a delta function peak at k — 

and k = (— if there was an ordering of the type 
suggested in On calculating the structure factor for 
various (ki,k2) at zero temperature, we find that S{k), 
apart from some fluctuation all through, has two clearly 
visible peaks at wave vectors (x'~T") ("T"' T" )' 
see Fig. [SJ However the height of these peaks are only 
about three times the average value, and they do not 
become sharper with system size. Thus, we find no ev- 



A 



V 




1 10 
r 

FIG. 7: (Color online) Graph showing the finite temperature 
correlation function ((/(O) — /(r))^) versus distance, r for var- 
ious values oi P — . The correlations are logarithmic at all 
temperatures, with the coefficient of log(r) varying between 

(2.45 ± 0.05) X 10"^ ~ III at /3 = and (4.12 ± 0.05) x 10"^ 
at /? = oo. 



idence of even incipient long-range order (hexatic-like, 
with power-law decay of the two-point correlation func- 
tion) in the system at T = 0. 

We also computed correlations of the / and (/) fields at 
various temperatures. At the end of each Monte Carlo 
step the 4> field was generated from the spin configura- 
tion by solving the discrete Poisson equation on the tri- 
angular lattice (Eq. (IA4[) ). This was done by invert- 
ing the Poisson equation in Fourier space as shown in 
Eq. (|A5p . The Fourier transforms were calculated us- 
ing the fast Fourier transform code provided in [26j . 
The spin configuration and the (j) field was then used 
to generate the /-field using Eq. p^ . In Fig. [Bl we 
have plotted the zero temperature correlation function 
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{{fLm — fi',m'Y) = ((/(O) — /(''))^) versus logr where r is eigenvalues A of Hi,(^i^„i) satisfy the relation 
the distance between the two sites. We see that this cor- 
relation function increases logarithmically with distance. 
Note that a logarithmic dependence of this function im- 
plies a 1 /r^ dependence of the bond-energy bond-energy 
correlation function (Eq. (|33|)). Fig. [7] shows the correla- 
tion function ((/(O) — /(r))^) at vaious values of /3. These 
correlations vary as log(r) at all temperatures with the 
coefficient varying between (2.45 ± 0.05) x 10^^ ~ 
at /3 = and (4.12 ± 0.05) X 10~2 at /3 = oo. 



V < (J5)^«(,„+i)^ + + s:(,„/). (40) 

This is true for all eigenvalues s^^i .^^^^ •Sa(i-i,m-Hi)' 
*a(i m) ^^"^ hence is valid as an operator inequality 



VII. THE LARGE-S QUANTUM KITAEV 
MODEL 



In this section we would like to discuss the ground- 
state energy of quantum spin-S* Kitaev model using a 
variational approach, which does not suffer from diver- 
gences. We obtain upper and lower bounds on the ground 
state energy. Our quantum mechanical Hamiltonian is 
normalised by the size of the spins. So 



^b{l,m.) - 

^2 



(JS) (S'a(;,m+1)) + iSa{l-l,m+l)) + ("^ a{l,7n) ) 



(41) 



(42) 



Therefore, for any wavefunction of all the 27V spins 
on the lattice we have 



(jsy {{s 

a 



(43) 



^ ' aeA 

Using the operator inequality AB > —{A"^ + B^)/2, 
where A and B are any commuting Hermitian operators, 
it is easily seen that H satisfies the lower bound 



E ground > "JN. (37) 

A better bound may be proved as follows. We write H 



H 



1 V 

S{S + 1) ^ ' 



(38) 



where Hm „i) is a 4-site spin Hamiltonian containing only 
the couplings of the site 6(Z, m) on the B-sublattice and 
its neighbours 



b(l.m) — J{Sa(l,m+l)^b(l,rn) + 



qy qy 

"^g(/-l,m+l)'^fc(;,m) 



^a(l,m)^b(l,m)) 



(39) 



The operators -fff,(/^m) can be diagonalized in a Hilbert 
space of 4-spins, i.e. a (25 + 1)^ dimensional Hilbert 
space. We note that S^(^i^„^+i)^ ^ii.m) are 

operators that belong to different sites, they commute 
amongst each other and with Hi,^i m) ■ Hence we can move 
to a basis in which these are diagonal. We now work 
in the subspace in which the basis vectors are eigenvec- 



tors of S' 



IX qy 

g(i,m+l)' *-'g(i-l,m+l) 



and S, 



with eigenval- 



ues Sa(i,™+i),Sa(i^i,m+i)'Sg(i,m) respectively Thus the 



where ((5^(,,„+i))') = {M%,m+i))^\^) and so on. Us- 
ing the fact that (lAI-ffK;,™) IV')^ < IV') and tak- 
ing the square root we get 



(44) 

This immediately gives 

(mm > E 

(;,m) 



(45) 

where the sum is over all the sites of the B-sublattice. 
Note that the terms in each of the squareroots are all 
real numbers. 

Using equation ((20| in the equation (|45p and observ- 
ing that for any site on the A-sublattice (<S'|^(;m))^ + 
+ iS:iUn) f = SiS + I) we get 



{ip\H\il;) > ~JN 



S 



S + 1 



For large S, 



(46) 



(47) 



is a lower bound, which shows an increase in the ground 
state energy due to quantum fluctuations. 

We now describe a variational upper bound for the 
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0.02 



o 
+ 




_^ 0.374/S — ■ 

0.02 0.03 0.04 
1/S 



FIG. 8: (Color online) Plot of E^ar + 0.5 as a function of ^ 
displaying a slope of 0.374 ± 0.005. 



ground state energy of the quantum model. 

We take a trial wave function which is a direct product 
of two-site pair wavef unctions: 



W = ^{l,m)\Tpail,m)b{L: 



(48) 



Each of these two wavefunctions \ipa{i.m)b{i,m)) is the 
ground state wavefunction of the two site Hamiltonian 
^2 given by 



— -[KSa{l.m) + ^b(l.rn)) + ^a(l,m)^b(l,m)]- (49) 

Using A as the variational parameter we calculate the 
resulting ground state energy of the full Hamiltonian. 

Let 7^ (A) = (5,^(,,„,))«., and e-(A) = 
(^a(i,m)^bii,m))'H2, where ( )'H2 denotes average un- 
der the ground state of [Eq. dH])]. Let E{X) be 
the minimum eigenvalue of Then, clearly for the 

wavefunction lib), we have 



(50) 



It is straightforward to determine the minimum eigen- 
value £'2 (A) by numerical diagonalization of the corre- 
sponding (25* -1-1)^ X (2S'-|-1)^ dimensional matrix for dif- 
ferent values of A. We can then also determine 7(A) and 
e^^(A) numerically from the corresponding eigenvector. 
The resulting upper bound on the ground state energy 
per site E^^^, was determined for different values of S by 
minimizing over A. Fig. © shows a plot of {Eg^ + 0.5) 
versus ^ . We see that the energy excess over the classical 
ground state energy varies as 1/5 for large S, 



Eg 



< E^, = --+ 0.374- + 0(-.) 



1 



1 



1 



JN 



g2' 



(51) 



The value of the coefficient of the 1/S correction term 
obtained by us (0.374) should be compared with the esti- 



mate 0.289 by Baskaran, et al. [B'l using the quadratic ap- 
proximation. The latter underestimates the true answer, 
as the modes that have frequency zero in the quadratic 
approximation would actually have a finite contribution. 
In our calculation, these corrections are not ignored. 
Also, the variational estimate can be improved with bet- 
ter choice of trial wavefunction [ip)^ say, in terms of six- 
site clusters. 

The ground state manifold of the classical problem 
has several soft modes. This suggests that one can find 
local ground states such that forming a product state 
with these as a basis yields a good approximation to the 
ground state of the full system. Then forming a linear su- 
perposition of such states locally, and forming products, 
lifts the degeneracy partially and gives a better varia- 
tional ground state. We perform such a calculation as 
follows. 

The sites of the hexagonal lattice can be divided into 
disjoint hexagons. Thus the full Kitaev model Hamilto- 
nian consists of interaction terms between sites within 
the same hexagon, and between different hexagons. We 
write H = J2hex Ethex + J2mt ^irit, whcrc the second sum 
is over the interconnecting bonds. Now some of the inter- 
connecting bonds are of X-type, some Y-type and some 
Z-type. We can redefine the directions x— , y— and z— 
at each site so that all the interconnecting bonds are of Z- 
type. Then the couplings between sites within a hexagon 
are alternating X- and Y-type, and the Hamiltonian Hhex 
for a single hexagon, with sites labelled 1, 2, ... 6 is of the 
form 



Hhe 



J 



S{S + 1) 

QX QX 
O5 Dq 



SlSf]. 



D2 O3 + D4 



ay qV 



(52) 



For the classical Hamiltonian Hhex, a state in which 
all spins are aligned in the same direction (j) in the xy— 
plane are classical ground states. The quantum state 
corresponding to this state denoted by \(j)4)4>4)(j)(f)) may be 
written as a product of single-site coherent states at the 
six sites 



5, (53) 



where represents the wavefunction of the spin at site 
J, polarized in the direction in the x-y plane, and can 
be written as 



(e'*S-) 



(54) 



where S = — lS^ is the angular momentum lowering 
operator at the site j, and \S)j is the state with — S. 

Since states for different (}> are classically degenerate, a 
better quantum-mechanical trial state is a linear super- 
position of these states 
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IV'i) = 



27r 



Ae^'^^-^l^^^). (55) 



The full wavefunction is a direct product of wavefunc- 
tions for different hexagons. The matrix elements of the 
spin operators are easy to calculate [1^ . In this state, the 
expectation value of Hint vanishes. Therefore we have 



1 {^l\Hhe.\i^i) 



6 (V'llVi) 

The expectation value of Hhex is easily seen to be 



(56) 



{lpl\Hhex\lpl 



'3J 



11 
US 



(57) 



which is a bit better than the variational bound for en- 

s 



ergy — SJ^j^, using the state 



VIII. DISCUSSION 

We have shown that the Kitaev model with classical 
spins shows no order-by-disorder, while there are plau- 
sible arguments that the quantum model does. This is 
because the mechanism of order -by- disorder in classical 
and quantum models is somewhat different. 

Consider a classical model, whose ground states form 
an M-dimensional manifold. Let G be one of the ground 
states. We expand the energy in the coordinates orthog- 
onal to the manifold, and look for small perturbations 
about the ground state. Keeping terms in the deviations 
from the ground state to quadratic order, and going into 
the normal mode coordinates, we get a quadratic approx- 
imation to the Hamiltonian in the transverse coordinates 



-m,{G)u,{Grq^l (58) 



where the sum over j extends over the n transverse de- 
grees of freedom. Then, the corresponding quantum me- 
chanical partition function in the quadratic approxima- 
tion is 



-'quad 



(59) 



For low temperatures, the G having the minimum 
value of the effective quantum mechanical free energy is 
obtained by minimizing the "zero-point energy" fq{G) = 
■^hujj. However, the classical partition function cor- 
responds to the case <^ 1, and at low temperature T 
is easily seen to be proportional to T"^/ Y[^j{G)- Thus 



the relative weights of different points G on the manifold 
are determined by an effective free energy fci{G) propor- 
tional to logWj. Clearly fq and fd are quite different, 
and states which are favored by one need not be favored 
by the second. In particular, fd depends more sensitively 
on the low frequency modes. 

If some of the oj^'s are zero, in this approximation, 
the classical partition function diverges, but the quantum 
weight has no singularity. This problem of zero modes 
also occurs in the calculation of [sj. In fact the zero 
frequency eigenvalue has a large degeneracy. 

More generally, finite h corrections in a quantum me- 
chanical system correspond to a finite temperature clas- 
sical model but in one higher dimension, and can be qual- 
itatively different. A simple example of this is a system 
of masses coupled by nearest neighbour springs in one 
dimension. In the classical case, the variance of displace- 
ments of masses at distance R varies as TR for large 
-R, and small T, but this quantity grows only as log R, 
both in the quantum case at zero temperature, and the 
classical case in two dimensions. 

In the path-integral formulation, the large-S quan- 
tum Kitaev model becomes a set of classical spins on 
a 2 -t- 1 dimensional lattice, with Kitaev couplings in 
two spatial directions, and ferromagnetic couplings in 
the time/inverse-temperature direction. It is quite plau- 
sible that in this 3-dimensional model, there is long- 
range order for low "effective temperature", but in the 
2-dimensional classical Kitaev model, the destablizing ef- 
fect of fluctuations is too strong. 

There are several interesting classical 2D systems, 
where thermal Order-by-Disorder is expected. The pro- 
totypical example is the system of Heisenberg spins on a 
kagome lattice, with nearest neighbour antiferromagnetic 
couplings. The expectation of order-by-disorder in this 
classical system comes from theoretical and Monte- Carlo 
studies, that suggest that at low ternperatures, the spins 
lie on a single plane as T — )■ [l^jfs^- This model can 
also be related to the height model at its critical point, 
within the quadratic approximation, suggesting that a 
single long-range ordered state (the a/S x ^/S state) is 
selected over the coplanar ones [3l|,[33. It would be in- 
teresting to identify the main reason for the difference in 
the behaviors in these models and the case studied here. 

We have mapped the finite and zero temperature prob- 
lem of classical spins onto a height-model interacting via 
an effective Hamiltonian Hsos- This effective Hamil- 
tonian depends on the temperature of the spin model. 
We have shown that the correlations of the height vari- 
ables vary as log(r), where r is the distance between the 
sites, for all temperatures of the spin model. The discrete 
height model with a pinning potential in 2D undergoes 
a roughening transition from a phase in which it is or- 
dered to one with logarithmic correlations between the 
height variables [11]. In the rough phase (T > Tr) of the 
height model, the coefficient of log(r) gives us a measure 
of the temperature [S^IjIIHI- So we see that the range of 
temperature [0, oo] of the spin model maps onto a range 
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of temperature of the height model which is in the rough 
phase. Note that as the temperature of the Kitaev Hamil- 
tonian is decreased, ((/(O) — /(r))^) in the SOS model 
increases. The increased fluctuations in / are accompa- 
nied by a decrease in fluctuations of the </> field, and the 
fluctuations of the spins {5} decrease with temperature, 
as expected. 



tential fields is 

4ia{l,m) + (j)a{l - 1, m + 1) + (t)a{l,m + 1) - 3(j)b{l, m), 

~ -~Qb(l,m) 

4>bil, m) + (pail + 1,™ - 1) + 4>a{l-, m - 1) - 3(/)a(/, m), 

= 
(A4) 
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where (f)a{l,m) = (j){a{l , m)) and (j)b{l,m) = (jj^alljUi) + 
Cz). Eq. (jA4[) can be inverted in Fourier space as 



-1 



9 - h{k)h*{k) 



3 h{k) 
h*(k) 3 





'Qb{k) 



(A5) 

where h{k) = 1 + exp[i(fc2 — fci)] + exp[ifc2]. Now, the 
difference in the / variables along the z axis is given by 

/(; -f i?,m+ 1) - /a,TO+ 1) = /fl - /o - 



Appendix A: / / Correlation functions at Infinite 
Temperature 



In this Appendix we calculate the asymptotic be- 
haviour of the / — / correlation functions at infinite tem- 
perature. There are two independent degrees of freedom 
at each A-site. We can choose these to be e(?,m;z) and 
e(Z, rn; y). Then e{l, m; x) is clearly —e{l, m; z) — e{l, m; y). 
At infinite temperature 



{e{l,m; z)e{l' ,m'; z)) = —Sij'Sm.m', 
45 

_2 

{e{l,m;y)e{l',m';z)) = — (5;.i"^m,m'- (Al) 
45 



i{l + r, m; z) + ^[0h(^ + r,m) ~ cj)a[l + r, m)]. 



(A6) 



We can write this in terms of the Fourier components as 
follows (taking (l,ra) = (0,0)): 



[fR- h] = — ^^[age(fc;z)-H/?je(fc;i/)]exp[-ifcir], 
h{k) - 3 



where — 1 + 



and /?g 



9-h{k)h*{k) 
h(k) - 3 



;i - exp[-jfc2]) 



— (expfifcil — l)(exp[— ifc2l). 
Q~h{k)h'{k) ^ 



(A7) 



For a field ra) on the lattice, we define the Fourier 
and inverse Fourier transforms as follows 

<t>{k) = — i=^exp[iA?.r](/i(f), 

f 

m = -i= exp[-zfc.r](/.(fc). (A2) 

k 

The vector r denotes the point (Z, m) in real space and 
k = {ki,k2) denotes the point (2^^' SttIvt) Fourier 
space, with k.f = ^ + The charge Qb(i,ni) at 

each B-sublattice site is given by 

Qb(Um) = -<^{l, m; z) - e{l - 1, rn + 1; y) 

+e{l,m + l;z) + €{l,m + l;y). (A3) 

The discrete Poisson equation that determines the po- 



Summing over r first and using Eq. (|Aip 



{{fR - fof)p=o = ^ 



l-exp[-ifci(i?-hl)] 



1 — exp[— ifci] 



[((a,,e(fc;z)-t-/3,,e(fc;y)) )]. 

(A8) 

The term within the square brackets can be shown to be 
equal to ^{1, r — i-gQ^^i — _). Thus the correla- 

^ 45 ^ o — cos /ci — cos ^2 — cos(A:i — A:2) 

tion function simplifies to 



{{}r ~ fof)^=o = 
r, n/T 2-^ 



1 -cos(fci(7?-|- 1)) 



45 LM ^ 3 — cos k\ — cos k^ — cos(fci — k-i 



(A9) 



We note that the factors of 1 — cosA:i cancel out in 
the expression. In the limit L, M 00 this summation 
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becomes an integral giving is equal to -^^^log[R] at large R [25j. Thus we have 

((f--f^)%-o-l,Gi^'^h (AlO) ((/,-/„n,^„.^log[i?]+0(l) for large R. 

457r 

where G(0,R) is the lattice Green's function on the tri- (AH) 
angular lattice between the points (0, 0) and (0, R) which 
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